#####################################
# Application
# Government formation
# Q2
#####################################
rm(list=ls())
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
source("cbq_function.R")

# load data
dat0 <- read.table("formation_new.tab",sep = "	",header=T)
dat <- dat0[order(dat0$case,dat0$realg),]
y <- dat$realg
Y <- ifelse(y==0,-1,1)
N <- length(y)
ind <- as.integer(as.factor(dat$case))
n_grp <- max(ind)
X <- as.matrix(subset(dat,select = c(minor,
                                     minwin,
                                     numpar,
                                     dompar,
                                     median,
                                     gdiv1,
                                     mgodiv1,
                                     prevpm,
                                     sq,
                                     mginvest,
                                     anmax2,
                                     pro,
                                     anti)))
n_coef <- dim(X)[2]
dat$ind <- ind

####################################
# estimate with CBQ
ms_cbq_stan <- cbq(N = N,
                   n_covariate = n_coef,
                   X = X,
                   Y = Y,
                   N_indx = n_grp,
                   ind = ind,
                   qtl = 2,
                   nchain = 5,
                   niter = 1000,
                   thin = 1)

save(ms_cbq_stan,file="cbq_q2.RData")
